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Summary 

For the recapture debugging design introduced by Nayak (1988) we consider the prob- 
lem of estimating the hitting rates of the faults remaining in a system. In the context of a 
conditional likelihood, moment estimators are derived and are shown to be asymptotically 
normal and fully efficient. Fixed sample properties of the moment estimators are com- 
pared, through simulation, with those of the conditional maximum likelihood estimators. 
Properties of the conditional model are investigated such as the asymptotic distribution of 
linear functions of the fault hitting frequencies and a representation of the full data vector 
in terms of a sequence of independent random vectors. It is assumed that the residual 
hitting rates follow a log linear rate model and that the testing process is truncated when 

♦This research was supported by the National Aeronautics and Space Administration, 
under contract NAG- 1-835. 


1 



the gaps between the detection of new errors exceed a fixed amount of time. 

Key words : Recapture Sampling; Software Reliability; Fault Hitting Rates; Asymptotic 
efficiency; Conditional Likelihood; Average Information; Interval Censored Sampling. 

l. Introduction 

Recapture sampling was introduced by Nayak (1988) as a way to get extra information 
for estimating the number of faults remaining in a system. By placing counters in the 
software we observe, in addition to the usual sequence of failure (i.e., error detection) times, 
the hitting frequencies of detected faults. Nayak’s (1988) discussion concerns the Jelinski- 
Moranda (1972) model and procedures useful for estimating the number of remaining 
faults. 

The present paper deals with the problem of estimating the residual fault hitting rates 
under a more flexible form of the model. Let A,- be the hitting rate of the faults remaining 
in a system when * — 1 faults have been corrected. Since eliminating faults improves 
reliability, we have Ai > A 2 >. . . and the effects of correcting the first, second, etc., 
detected faults can be defined by £, = A, — A, + j. To model the failure times Ti,T 2 , ..., we 
assume that the failure gaps Y, = T t - T,_i, (T 0 = 0) are independent and exponentially 
distributed with rate parameters A,-; note that Y{ can be interpreted as the smallest of the 
times of encountering the remaining faults. As for Nayak’s (1988) model, counts {M,(t)} of 
repeated error occurrences are assumed to be independent homogeneous Poisson processes 
with rate parameters 

In attempting to get consistent estimators of the A,-, the following difficulties arise: 
(t) the likelihood function based on the ordered failure times and hitting frequencies is not 
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indexed by a measure of the amount of information and (it) the number of parameters 
increases in proportion to the amount of information. 

In view of (it) we study estimation under the following model: 

A i = , a, 0 > 0 (l) 

= a«“* (<_1) (l - e~ p ) i — 1,2, ... 

This model yields a tractable form of the likelihood function and is considerably more 
flexible than Jelinski-Moranda (1972) model. The latter model A, = {u — i + 1 )</>,i = 
l,2,...,t/, 4> > 0 assumes a common rate, & = <f>. Empirical evidence indicates, however, 
that faults may have different hitting rates (Nagel, Scholz, Skrivan, 1984). The effects of 
correcting faults will typically decrease since faults having the highest hitting rates are 
likely to be detected early during the debugging process. 

The full likelihood function is presented in Section 2. This likelihood, however, is 
not indexed by a measure of the amount of information so we consider in Section 3, a 
conditional likelihood obtained by conditioning on the observed number of detected faults. 
In this case the sufficient statistic can be represented as a sum of independent noniden- 
tically distributed random vectors. Using this representation, we show in Section 4 that 
certain moment estimators are asymptotically normal and fully efficient. In Section 4 we 
also consider estimating on upper bound for A r +i where r is the number of faults even- 
tually corrected. Fixed sample properties of the moment estimators and the conditional 
maximum likelihood estimators are evaluated, through simulation, in Section 5. 
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2. Notation and other preliminaries 


Suppose a system is tested until the failure gaps first exceed a fixed amount of time s. 
We observe T U T 2 , ...,Tr where R , the number of faults eventually detected, is determined 
by truncating T, = 7\ - T,_i , *'=1,2,... over the interval (0, s). That is, data is obtained 
through an interval censoring model in which we observe the realization R = r providing 
Yi < s,i = 1,2, ..., r and Y r+i > s . The total test time r = ]T)f Y i + 5 is then random 
while for Nayak’s (1988) censoring model, r is a fixed time of testing. 

Since yi,y 2 >--> are assumed to be independent exponential random variables, the 
joint density function of (R, Yi,Y 2 ,. . . , Yr) is 

r 

/(yi,y 2 , = exp(-A r+ is) A t exp(-A t y<), 0 < y, < s,i = 1,2, . . . , r (2) 

t=i 

Explicit formulas for the case r = 0, will be omitted since if r = 0, there is little information 
available for inference. 

Recall that counts {M, (f)} of repeated error occurrences are assumed to be indepen- 
dent Poisson processes with rate parameters £,. For the fault detected at time f, let M, 
be the n umb er of times this fault is accessed during the interval Given t\,t 2i ..., t r 

the vector (Mj, Af 2 , ..., M r ) consists of independent Poisson random variables with means 
= 1, 2, ...,r. 

The joint density of (R,Y\, Y 2 , —,Yr, Mi,M 2 , ...,Mr) is 

r r r r 

exp(— A r+ is) A, exp( — A,y t ) JJ[&( X* + s )] m ‘ ex P[ - f*( X* yy+ 5 )]/ m ‘- 

t=i »=i y=*+i j=t+i 

= exp(-Air) JJ A,-[£,( X] yy + s)] m 7 m t ! ( 3 ) 

*=1 j=i + 1 
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Under (l) the likelihood function given by (3) is maximized by 


a= (r + ^m,)/r 
1 

r r 

$ = ln{ 1 + ^2 m */ “ 1 ) m * + r ( r “ 1 )/ 2 H' 

l l 

However, since this likelihood is not indexed by a measure of the amount of information, 
the asymptotic properties of these estimators cannot be deduced from the form of the 
likelihood. 


3. A Conditional Likelihood 

Under the censoring model of Section 2, the conditional density of Yj, Y 2 , • ■ • , Y r , given 
R = r, is 

r 

1 1 exp( — A,y t )[l ®xp( ^t®)] » 0 < y t ‘ < s; i — 1, 2, . . . , r (4) 

t=i 

Consider the sequence defined by 

Z 1 = (Y l ), Z fc = (y fc ,Af lfe ,M 2fc ,...,M (fc _ 1)fc ) (k = 2,3, ...,r + l) (5) 

where M{k,i < k, is the number of repeated occurrences of error t during the interval 
(Tjt_i,Tfc] and, for convenience, we let Yi-+i = s. The last interval ( T r ,r] has fixed length 
s while (Tfc_i,Tjt], k < r, has random length Y k . Our earlier assumption that are 

independent homogeneous Poisson processes implies that {M^}, t < k,k = 2,3 ,...,r + 1 
are conditionally, given Y\, Y*j, ...,Yr andiZ, independent Poisson random variables with 
means £ t y*. 
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Since counts of events occurring in different intervals are independent and the interval 


lengths are also independent, it follows that Z\, Z 2 , Z r + \ are conditionally, given R — r, 
independent with densities 

9k{Vk\ ™ lk ,m 2k , • • • , ra(fc-i)fc) 

k-i 

= X k e- x ^{l-e- Xke )- 1 Y[{tm) mik e-i iyk /m ik \ (k = 1,2, ...,r) 

1=1 

r 

= n^‘ 5 ) m * <r+i> e_€,a / m *( r+i ) ! ( fc=r + i ) ( 6 ) 

t=i 

Substituting our model for A< and & in (6) and simplifying by using £ l=1 & = 

gives the conditional log likelihood 

fc-l k-l k- 1 

Ik = C k (a,0 ) + /n(l - e - ^) m tfc + In a ^ m,* - /?^(t - l)m«fc - ay* + C (7) 

t=i t=i »=i 

Here C does not depend upon (a,/?), and 

C k (a,/3) — Ina — (3{k — 1) — ln{ 1 — exp[— ase - ^* -1 )]} (k = 1,2, ...,r) 

= ctse~P r (k = r + 1) 

In terms of A*, £*(«,/?) = /n[A*(l - e -A * 4 ) -1 ], k = 1,2 and C r+ i(a,/3 ) = A r+ is. 

Let Vjfe = (V r u,V 2fc ,V 3 *), fc = l,2,...,r + l be defined by 

fc-i fc-i 

V lk = Y k , V 2k = '£M ik , V 3k = ^2(i-l)M ik ( 8 ) 

»'=i i=i 

where, as before, Y r+i = s. In reference to l k in (7) and the full conditional likelihood 
l c = l k note the following: 

(i) Vj , V 2 , Vr+i are independent. 

(ii) Vjt is a minimal sufficient statistic for l k . 
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(iii) S r = Z[ +l Vk is a sufficient statistic for the family defined by l c . 

(iv) S r = (5i r , S 2r , S 3r ) is given by 5 lr = r, S 2r = S 3r = Z)i(* “ 


To study the asymptotic distribution of S r , we shall need the following : 

THEOREM 1. Let Vk be defined as in (8) where (Yk,Mik,M 2 k, —, M(k-i)k), k = 1,2, 
. . .,r + 1 are independent with densities (7 ). Let the means, variances and covariances 
of the elements of V* be denoted by pik = E(Vk) ] * = 1,2,3 and Oijk = Cov(V{k,Vjk); 


i,j =1,2,3. Then 

Mi* = A* 1 - s(e Aka - 1) _1 

M3* = Ml *<** 

<^12* = ^llk^l — ^*) 

^22* — Ml*(^l — ^*) + <7ll*(Ai — ^*) 2 

<^33* = Ml *&* + ^11*°* 

*- 1 

a* = B*- 1 * 

t=i 

Ml(r+1) = 5 


M2* - Ml*(Al - At) 

<r„t = Aj 2 - S V>*(e i ‘*-l )- 2 

<713* — <^11*0* 

<723* = Mi * a * + <7n)fe(Ai — Afc)a^ 
bk = £(*-D 2 * 

i=i 

<7ll(r+l) = 0 


PROOF: The mean pik and variance 0nk of Y* can be obtained directly from the fact 
that Yk has an exponential distribution truncated over the interval (0,s). Since V 2 k is 
conditionally, given Yk, a Poisson random variable with mean Yk J2i=i £* = ^*(^1 — ^*)> 
and Var(V 2 k) = E[V ar(P 2 *|i*)] + Var[E(F 2 fc|Ylk)] we immediately have the expression 
given for o 22 k- Similar calculations can be made for by noting that is a linear sum 
of the conditionally (given Yk) independent Poisson random variables < k. These 

being routine calculations, further details are omitted. 
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The Average Information Matrix: 


The second and mixed partial derivatives of /* are as follows 
= -(1 + E + 5 J £>- 2 A Je A »'(e A *' - 1)“ 2 


^5 = «*-'(* - l)A t (e A ‘* - I)' 1 - s 2 a~ 1 (k - l)A! e A *V*' - 1)“ 2 
oaop 

~ = -s(k - l) 2 A fc (e A ** - l)-‘ + s 2 (k - l) 2 Aje A * V - l )’ 2 - e V ~ l )" 2 E 

°P i=l 

The Fisher Information Matrix A* = (a t y*), based on /*, is given by 

o2 / 

anfc = — £( g- = a Vifc + a 2 (A*<7n* — A*/ii*) 


°12fc — ~ ~ ^) a 1 (Afc <J llfc Afc/ijfc) 


a 22 * = = (^ - !) 2 (AfcCTufc - Xkfiik) + (Ai - Afc)^ifce^(e /? - 1) 2 

To simplify the calculation of these quantities, write the second derivative expressions 
in terms of the moments of Y* . For example, 


s 2 e A*«( e A fc « _ 2 = X k 2 — <7nifc 

0 = -a-*‘£ M ik -a~ 2 Xla nt 

t=l 

Then am is obtained by using first the moment of V 2 * given in Theorem 1. Similar 
calculations yield the expressions given for ai 2 * and a 22 * • 

Since the distribution of Y* converges (as k — * oo ) to a uniform distribution on the 


interval (0, s) the moments of Y* converge to the corresponding moments of the limit- 
ing uniform distribution (Serfling, 1980 p.14). We, therefore, have lim* —*■00 Ml k — 
limjfc^ootfiiJfc = s 2 /12, and lim*^ ( k - l^A* = 0 for 7 = 0,1,2. Using these limits to 
get the limiting average information matrix A = lim(l/r)(Ai + A 2 + • . . + A r+ i), we have 
A = (aij) where an = s/2a, ai 2 — 0, a 22 = ase^(e^ — l) -2 /2 
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4 . ESTIMATION 


Let a and (3 denote estimators of a and /? defined by 

r 

a — ^ Mi / r 

t=i 

r r 

1= 1 t=l 

Here we derive a and /? by equating (1/r) £| Mi and (l/r) £j(* — 1)-Wt to the appropriate 
mean vector elements of the asymptotic distribution of (l/r)S r . 

THEOREM 2 . Under the assumptions of Theorem 1 , (1 /r)S r has a limiting (as r tends 
to infinity) normal distribution with mean vector p' = (pi , /* 2 > ^ 3 ) and covariance matrix 
(l/r)E. £ = where 

Pi - s/2 P2 = as/2 pz = ocs(e@ — l ) -1 /2 
on = s 2 / 12 oi2=as 2 /\2 013 = as 2 (e^ — l ) _1 /12 

<722 = <*>s/2 + a 2 s 2 / 12 

<723 = — l) - 1 [s /2 + as 2 / 12 ] 

<733 = (<m/ 2 )[(«* - I )" 1 + 2 (e^ - I)" 2 ] + (a 2 s 2 / 12) (e* - 1)~ 2 

PROOF: The elements of p and E are given by pi = lim(l/r ) YLk=\ = 1» 2 ,3 and 

Oij = lim (1/r) £fcti °ijk as r tends to infinity, where the terms in these sums are the 
moments given in Theorem 1. Since and <7,;jt converge to finite limits as k tends to 
infinity, we have pi = lim Pik and <7 t y = lim Oijk- Thus the calculations are similar 
to those discussed at the end of Section 3. The remainder of the proof requires showing 
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(Serfling, 1980 , p. 30 ) that 


r-f -1 

lim(l/r) ^2 h kr = 0 (9) 

k=\ 

where 

h kr = E\U k I(U k > e 2 r)} 

3 

Uk = Y,( Vik ~ fii tf (10) 

t=i 

and /(.) is the indicator function. Since h kr < {t 2 r)~ l E{U k ), the limit in ( 9 ) can be 
established by examining the fourth central and product moments of the V ik . Further 
details are given in the Appendix. 

Note that a = yi(/ii,M2,M3) and (3 — <72(7*1,7*2, M3) where <71(21,22,23) = 22/21, 
^2(21,22,23) = ln( 1 + 22/23), and p' is the mean vector of the asymptotic distribution of 
(l/r)S r . Applying the 6 -method gives the following. 

COROLLARY 1 . Under the assumptions of Theorem 1 , a and (3 have a limiting (r — ► 00) 
normal distribution with mean vector (a, ( 3 ) and are asymptotically independent with 
variances o s, 2 = 2 a/rs, o^ 2 = 2(e P — 1 ) 2 e~^ jras. 

COROLLARY 2 . Under the assumptions of Theorem 1 , A g +i = ae~P q , for fixed q < r + 1 , 
has a limiting (r — ► 00) normal distribution with mean A ?+1 and variance 
o 2 = { 2 a/rs)e~ 2 ^[l + 2 <7 V " l) 2 ^]- 

By comparing the form of (a,/?) with that of the maximum likelihood estimators given 
in Section 2 , it is evident that the latter are not consistent. Consistency, however, depends 
upon the choice of index for the likelihood function. That is, there may be other ways of 
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conditioning (i.e., other families of conditional likelihoods) which imply the consistency of 
a and (3 . 

Corollary 2, gives an asymptotic basis for estimating an upper bound on A r+ i but 
does not justify directly estimating A r+1 . The hitting frequencies at the upper end of the 
vector (Mi, M 2 ,..., M r ) will tend to be small (often zero). In such cases, for a suitably 
chosen index q < r + 1, \ q+ i - X r +i = Yli= q +i £* wil1 tend to be sma11 so that ^<?+i ma y 
be a close upper bound on A r+ j. 

5. Efficiency and Bias 

Since a 5 ? = (ra u ) -1 and a ^ = (ra 2 2) -1 where a u and a 22 are given at the end of 
Section 3, it follows that a and f3 are asymptotically fully efficient. 

To study the fixed sample properties of a and /?, we simulated their values under 
the conditional model defined in (6). This was done by generating 1,000 replicates of 
Z = (Z 1 ,Z 2 ,...,Z. +1 ) for the values of r shown in Table 1. In addition the conditional 
maximum likelihood estimates a c and /3 C were calculated for each replicate by maximizing 
l c = X!i +1 Ik where Ik is defined in (7) . 

Table 1 shows the bias and mean square error (MSE) for a and /? and also the bias 
and MSE for a c and j3 c . The conditional MLE clearly have smaller bias than a and (3 and 
thus would probably be preferred in applications. 

6 . Final Remarks 

We have presented procedures that may be useful for estimating an upper bound on 
the hitting rate of the faults remaining in the system. There does not seem to be any 
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published data where software testing counters have been used, although Nayak (1988) 
credits this testing method to Huang (1977). A few studies have used replications of a 
debugging sequence and multiple versions of a program to automate the process of error 
detection (e.g., Nagel, Scholz, and Skrivan, 1984). 

The model in (1) and, more generally, any model in which the £, are defined by 
£ = A, — A, + i, assumes that the effects of correcting faults are additive. If for some input 
condition, two or more faults give an error in the output, then the additivity assump- 
tion does not hold. Although disjoint fault sets are more common, cases of a nonempty 
intersection have been observed experimentally (Nagel, Scholz, and Skrivan, 1984). The 
additivity assumption can be preserved by (o) labeling such faults as a single fault and 
recording only the hitting frequency of the union of the fault sets or (6) partitioning the 
union of these sets into disjoint regions and recording the hitting frequencies as if each 
region corresponds to a different fault. The former is a simpler and more practical method 
since the hitting frequencies of the common fault regions may tend to be small. 

In analyzing data obtained by recapture sampling, it is important to use estimators 
that are consistent (i.e., converge in probability to the true parameters). Although consis- 
tency is not guaranteed for estimators based on an unconditional likelihood, we have shown 
that certain moment estimators, obtained in the context of a conditional likelihood, are 
consistent and asymptotically normal; it is likely that the conditional maximum likelihood 
estimators are also consistent. Our asymptotic evaluation of the estimators assumes that 
the number of detected faults is large. It is, therefore, of interest to determine whether 
consistency of these or other estimators will hold when conditioning on other quantities, 
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such as the total number (tto + r where m — ^») of error occurrences. 

7. APPENDIX 

Before we give the proof of Theorem 2, we require the following moments, which 
can be obtained by using (8) where V 2 i t and Vjjjt are linear function of the conditionally 
independent Poisson random variables {M,*}. 

E[{V 2 k ~ M2fc) 4 ] = 3 alE{Y?) + a 0 Mife + 6 alE[Y k {Y k - fi lk ) 2 } + a A E\{Y k - Mifc) 4 ] 
E[(V 3^ — fi 3 k) 4 } = 3 a 2 E(Y£) + a^n\k + 6 a,iE[Y k (Y k — Mifc) 2 ] + a i E[(Y k - Mifc) 4 ] 

— A t 2fc) 2 (^3fc — M3fc) 2 ] = <A a o£tik + <t> a oE[{Y k — Mifc) 3 ] + <£ a oMifc 

+ <f>a^ lk E[{Y k - /x lfe ) 2 ] + (a\/al)E[(V 2k - » 2k ) 4 } 

k - 1 

«P = E(*-1) P & p = O’ 1 ’ 2 ’ 3 ’ 4 * = a o 1(12 - “o 1 2 °1 

t=l 

The limits, as k — ► oo, of these moments are finite because (*) the Y k are bounded and 
converge in distribution to a uniform distribution on the interval (0, s), and («) for each 
p=0,l,2,3,4, a v = a(l - e _/? ) ^“/(i - converges to a finite positive limit as 

k — y oo. Therefore, E(U%) converges to a finite limit where U k is defined in (10). Thus 
as r tends to infinity, the limit in (9) is zero, which proves Theorem 1. 
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r 


A. 

OC c 

A 

Pc 

a 

P 

15 

Bias 

0.000066 

0.005990 

-0.011208 

0.055903 


MSE 

0.000196 

0.003188 

0.000295 

0.006629 

20 

Bias 

-0.000001 

-0.000001 

-0.007733 

0.041084 


MSE 

0.000207 

0.002241 

0.000276 

0.004115 

25 

Bias 

0.000057 

0.000002 

-0.004581 

0.032183 


MSE 

0.000103 

0.000680 

0.000106 

0.001701 

30 

Bias 

0.000000 

0.000000 

-0.002472 

0.022305 


MSE 

0.000074 

0.000789 

0.000088 

0.001288 


Table 1. Bias and mean square error (MSE) of the conditional maximum likelihood esti- 
mators a c , f3 c and the moment estimators a, P based on 1,000 simulations with a = 0.10 
and (3 = 0.35. 


